%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lorentz model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
figure(1);          clf;

gam=2*pi*   6e6;
del=[-5:.1:5]*gam;
w0 =2*pi*500e12;
w  =w0+del;

% polarizability 

alp=1./(1i+2*del/gam);
amp=abs(alp);
pha=atan(imag(alp)./real(alp))/pi;

subplot(221),plot(del/gam,amp,'b',del/gam,pha+(pha<0),'r','LineWidth',1.5)
set(gca,'FontSize',16,'LineWidth',1.5)
xlabel('\fontsize{16}(\omega-\omega_0)/\gamma')
ylabel('\fontsize{16}{\itpolarizabilidade}')


% refraction index

nrf=1e20* (w0^2-w.^2)./((w.^2-w0^2).^2+gam^2*w.^2);
alf=1e20* 2/6e15 .* w.^2.*gam./((w.^2-w0^2).^2+gam^2*w.^2);
[mxm,ixd]=max(nrf)
[mnm,imd]=min(nrf)

subplot(222),plot(del/gam,nrf,'b',del/gam,alf,'r','LineWidth',1.5)
hold on;
subplot(222),plot(del(imd)/gam*[1 1],[-.001 mnm],'--k',del(ixd)/gam*[1 1],[-.001 mxm],'--k','LineWidth',1)
set(gca,'FontSize',16,'LineWidth',1.5)
xlabel('\fontsize{16}(\omega-\omega_0)/\gamma')
ylabel('\fontsize{16}{\itn-1}    ,   {\it\alpha}')
text(del(imd)/gam*1.5,-.00085,'\fontsize{14}\omega_2')
text(del(ixd)/gam*2.7,-.00085,'\fontsize{14}\omega_1')